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Abstract 

We demonstrate the universality of the spectral correlation functions of a QCD inspired ran- 
dom matrix model that consists of a random part having the chiral structure of the QCD Dirac 
operator and a deterministic part which describes a schematic temperature dependence. We 
calculate the correlation functions analytically using the technique of Itzykson-Zuber integrals 
for arbitrary complex super-matrices. An alternative exact calculation for arbitrary matrix size 
is given for the special case of zero temperature, and we reproduce the well-known Laguerre 
kernel. At finite temperature, the microscopic limit of the correlation functions are calculated 
in the saddle point approximation. The main result of this paper is that the microscopic uni- 
versality of correlation functions is maintained even though unitary invariance is broken by the 
addition of a deterministic matrix to the ensemble. 



1 Introduction 



The study of level correlations in quantum physics has a long history. It was first realized 
in nuclear physics that the spacing distribution of compound nuclear resonances is given by 
random matrix theory Q. Much later it was found numerically that the spectral correlations of 
systems as simple as certain quantum billiards |2|, [| or polynomial potentials with two degrees 
of freedom 0] can also be described by random matrix theory. The expectation that the spectral 
correlations of the invariant random matrix ensembles are universal is most clearly expressed by 
the Bohigas conjecture ||: The spectral correlations of a quantum system are given by random 
matrix theory if the corresponding classical system is chaotic. Although subsequent work has 
provided considerable analytical insight regarding this conjecture || ||, [7|, ||, ||], a proof has 
remained elusive. Meanwhile, the universality of random matrix correlations has been tested in 



a wide variety of systems including the zeros of the Riemann zeta function [10|, sound modes in 



crystals [11], and the eigenvalues of the QCD Dirac operator 12]. 

We address the question of universality within the context of random matrix theory. An 
observable, typically a spectral correlation function, will be called universal if it is stable against 
deformations of the probability distribution. Recently, this topic has attracted a great deal of 
attention. Two types of deformations have been considered: deformations which maintain the 
invariance of the random matrix ensembles [13, 14, 15] and those which violate that invariance 



i~G| , 17, 15, IS], We will consider spectral correlation functions for an example of the latter class. 
In all cases which have been studied, spectral correlations measured in units of the average level 
spacing are found to be universal even for deformations that change the spectral density on a 
macroscopic scale. 

Specifically, we investigate the stability of the spectral correlations of the chiral ensembles 
|20[| . They are constructed to describe the fluctuations of the Dirac eigenvalues in lattice QCD 
JI^ ] and are relevant to the theory of universal conductance fluctuations. Because of an under- 
lying Ua{^) symmetry (i.e. because the Dirac operator commutes with 75) , the eigenvalues of 
chiral matrices occur in pairs ±A. Therefore, one can consider three types of universal behavior 
each of which is given by invariant random matrix ensembles: (i) correlations in the bulk of the 
spectrum pl| , (ii) correlations near the edge of the spectrum |^], and (iii) correlations near 
A = 0. In this case, the microscopic spectral density (i.e., the spectral density near zero on the 



scale of a typical eigenvalue spacing) is universal. This has been illustrated for both invariant 



23] , 24] and non-invariant deformations [25] 



The invariant random matrix ensembles and the chiral random matrix ensembles are part 
of a larger classification scheme: Altland and Zirnbauer have shown that there is a one to one 
correspondence between random matrix ensembles and symmetric spaces [p6| . 

The microscopic spectral density is of immediate physical interest. According to the Banks- 
Casher formula pTjl , the spectral density at zero is directly proportional to the order parameter 
of the chiral phase transition in QCD (i.e., the chiral condensate). The microscopic spectral 
density provides information regarding the approach to the thermodynamic limit. This has 
been demonstrated, for example, in connection with the dependence of the chiral condensate 
on the valence quark mass pq ]. Recently, it has been verified by direct lattice calculations that 
the microscopic spectral density of lattice QCD is given by one of the chiral ensembles [29]. 



Moreover, the microscopic spectral density enter in sum-rules |30| for the inverse eigenvalues of 
the Dirac operator. 

The model which we shall consider is the chiral unitary ensemble perturbed by the lowest 



Matsubara frequencies, ±7rT. This model was introduced in [31] as a model for the chiral phase 



transition. Indeed, the average spectral density for this model undergoes a transition from one 



semicircle at zero temperature to two disjunct semicircles at high temperature. In [25] it was 
shown that the spectral density of this model on the scale of individual average level spacings 
is independent of the temperature. 

The aim of the present paper is to show that, below to the critical temperature, all spectral 
correlations measured in units of the average eigenvalue spacing are independent of the tem- 
perature. This will be achieved using the methods of super-symmetric integration introduced 
by Guhr (3^] and certain super-symmetric Itzykson-Zuber integrals 35, |36|, |37|] . At zero 
temperature our results coincide with earlier work pO, 38, 39, 24). 



The outline of the paper is as follows. In section 2, we introduce the random matrix model 
and express the correlation functions in terms of a partition function. In section 3, we reduce this 
partition function to an integral over a supermatrix, a, of much smaller size. The choice of an 
explicit parametrization for a is essential for the rest of the paper. Thus, we will argue in section 
4 that the correlation functions can be obtained by deforming the non-compact parametrization 
required for uniform convergence into a compact parametrization of a. This will be done through 



a detailed investigation of the one-point function. In section 5, we show that all correlation 
functions follow from a two-point kernel. This will be evaluated at zero temperature in section 
6, and we will reproduce the well-known Laguerre kernel. In section 7, we evaluate the two-point 
kernel at finite temperature using a saddle point approximation. Our primary result will be that 
the correlation functions are independent of temperature except for a trivial rescaling of their 
arguments. 

2 The random matrix model and the partition function for cor- 
relations 

The partition function of Euclidean QCD is given by 

Z E = {det(i 1 -D + iM))s E (1) 

where j-D is the Dirac operator, M is the mass matrix, and the brackets denote the average over 
gauge field configurations with respect to the Euclidean Yang-Mills action Se- Using a lattice 
regularization with anti-periodic boundary conditions in the time direction, we can expand the 
fermion fields as a sum over Matsubara frequencies to obtain 

N n 

^(x,r) = ^ ^(x)exp(^(2/-l)Tr) . (2) 

fc=l Z=-n+l 

Here, T is the temperature (i.e., the inverse length of the time axis), ./V is to be identified with 
the volume of space, and n denotes the total number of Matsubara frequencies retained in the 
expansion. The 4>k are properly normalized spinors. 

Using this basis, we separate the time derivative in the Dirac operator, 7o<9 T , from the 
remaining terms. In terms of the matrix elements of the Dirac operator the partition function 
for Nf massless flavors can then be written as 

Z 'gauge = ^det^ 

The identity matrix of size N is denoted by 1 n > and O contains the Matsubara frequencies B = 
diag (— (2n — 1)ttT, . . . , — ttT, ttT, . . . , (2n — 1)ttT). In (||), we have used a chiral representation 
of the Dirac matrices. For fundamental fermions with three colors, the matrix C is a complex 
matrix of size 2nN x 2nN. The detailed form of C depends on the particular gauge field 
configuration. 
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(3) 



The corresponding random matrix model is obtained by replacing the matrix elements of C 
by independent random variables with a Gaussian distribution. Thus, instead of (||), we study 
the properties of 



Z RM = N D jd[C]e- N ^ c1c det N f 



iCt \ / iQ® ljv 
iC ) [ i&®l N 



(4) 



where the measure d[C] is the Haar measure defined by \[ m n d(Re C mn )d(lm C mn ), and Njj is 
a normalization constant to be specified later. We wish to stress that this partition function 
is a schematic model of the QCD partition function. For example, we have ignored all spatial 
dependence of the matrix elements of the Dirac operator, and the critical exponents of this 
model are necessarily those of mean field theory. It is our claim, however, that this partition 
function belongs to the same universality class as the QCD partition function with respect to 
local spectral fluctuations. 

In this paper, we will study the spectral correlation functions of the model given by (||). For 
simplicity, we will restrict our attention to the case in which we retain only the lowest Matsubara 
frequencies, ±7rT. This enables us to replace by irT using a unitary transformation. After 
redefining N, the matrix C can be taken to be a complex N x N matrix. 



As was shown in [31, 25], this model shows a second-order phase transition at 7rTS = 1. 
Below this temperature, chiral symmetry is broken spontaneously with the chiral condensate 
given by 3 = — 7r 2 T 2 £ 2 ). At all temperatures, the spectral density follows from the 

solution of a cubic equation |3l], ^5|. In particular, we mention that, according to the 
Banks-Casher formula ]p7|] , the spectral density at zero virtuality is simply related to the chiral 
condensate: 

N ■ W 
In [25] it was shown that the microscopic spectral density, defined as 



is independent of the temperature parameter in the random matrix model. This result greatly 
adds to our understanding of the empirically observed universality of the microscopic spectral 



density in lattice QCD |28|, |2S[ and instanton liquid simulations [20]. Because of these results 



and numerical results for higher-order correlation functions in the bulk of lattice QCD Dirac 



spectra [O], we expect that the higher-order microscopic correlation functions are universal as 
well. 

The spectral density of a matrix with eigenvalues X k is given by 

p(x)=J2S(x-X k ) . (7) 

k 

The fc-level correlation functions are then defined as 

Rk(xi, ...,x k ) = {p(xi) . . . p{x k )) , (8) 

where (• • •) denotes the ensemble average with respect to (Q). However, in our approach is is 
more convenient to work with correlation functions of the resolvent 

where D is the random matrix Dirac operator 

and e is a positive infinitesimal. The spectral density is then given by 

p(x) = -ImG(x). (11) 

7T 

The corresponding correlation functions are defined by 

R k (x 1 ,...,x k ) = (l[-Ti 7 1 n \ • (12) 

The spectral correlation functions follow by taking the imaginary part of each trace. In the 
following, we will drop the ie\ term, assuming that xi has a positive infinitesimal imaginary 
part. 

Using the identity 

TrA- 1 = [det(A-j)/det(A + j)]\ j=0 , (13) 

one can express the correlation functions in terms of a generating function 



1 A d 



R k { Xl ,...,x k ) = J[—Z k (j 1 ,...,j k ) 



(14) 



jl=0 



Z k (3 1 ,... t3k ) = N D Jd[C\e det /fl g ___ . (15) 



The prefactor, Nd, in (|l5| ) is to be chosen such that Z k (0) = 1. In the rest of the paper, we will 
restrict our attention to the quenched problem with Nf = 0. This is due to the fact that the 
current method does not allow us to find a solution for arbitrary Nf. The reason for this will 
be discussed at the end of section 5. 

3 Reduction of the partition function 

The expression ( |l5| ) for the generating function of correlators involves an integral over a matrix, 
C, with 2N 2 de grees of freedom. When the deterministic part, 0, is absent, we can exploit the 
unitary invariance 

C^UCV~\ (16) 

with U and V unitary matrices, in order to rewrite the partition function as an integral over 
the eigenvalues of C only. It can then be evaluated easily using, for example, the orthogonal 



polynomial method [41, |3g]. The presence of O destroys unitary invariance, and standard 



methods are unsuitable. However, unitary invariance is not essential in the supersymmetric 
formulation of random matrix theory [42, 43]. Below, we adapt an approach based on the 



supersymmetric method which is described in [32]. This method exploits the determinantal 
structure of the correlation functions and allows us to calculate all /c-point correlation functions 
at the same time. Readers not familiar with the use of supersymmetry methods in random 
matrix theory are referred to [42, 43, 32]. 



In this section, we will express the partition function, (|T^), as an integral over a 4fc x 4/c 
supermatrix in which N appears only as an overall factor in the action. This is particularly 
useful for the investigation of the thermodynamic limit (i.e., N — > oo). 

The ensemble average in ( p~5|) is performed by writing the determinants as Gaussian integrals 



over commuting and anticommuting variables [43]. Equivalently, the product of determinants 



and inverse determinants in (|T^) can be considered as a superdeterminant and can be written as 
a Gaussian integral over a complex supervector. In order to do this we introduce the following 
supermatrices 

C = diag(C,..., C;C,..., C) , 
x = diag(xi, . . . ,x k ;x 1 , . . . ,x k ) , 

j = diag(-ji,...,-j fc ;j'i,...,j fc ) . (17) 



Here, a semicolon separates the k boson-boson blocks from the k fermion-fermion blocks. Each 
block is of size N x N . Further, x p and j p are understood to be multiplied by the N x N identity 
matrix. 

In order to exploit the chiral structure of the problem, we introduce a pair of complex 
supervectors, 4>a and 4>b, each of size (kN;kN). The two numbers refer to the number of 



written as 

k 



commuting and anticommuting components. Then, the product of determinants in (15) can be 
3n as 

det(D - x, + j { ) 

= su< 

d[(f) A }d[(j)B\ exp 



n<let(D-xi-ji) 



sdet 



-i 



-x + j C^ + ttT 
C + irT -x + j 



4 



-x+j Ct + ^T 
C + itT -x+j 



4>A 
4>B 



(18) 



where we have taken = irT. In ([Tg|), the blocks inside the superdeterminant refer to the chiral 
structure, i.e., each block is itself a supermatrix. The measure on the right hand side of (|l8|) is 



k 

n 

P =i 



.1=1 



2tt 



k 

n 

9=1 



■ N 

n 

.m=l 



gray 



(19) 



Here, <^ refers to either <j>^ or 0b; 0^ and (f)] jm denote the commuting and anticommuting 
components of 4>, respectively. The ranges of the indices follow from the limits in the products. 
The integrals with respect to the anticommuting variables are normalized according to the 
convention |f44}| , J d\X = 1- The constants in (|Tg|) are chosen so that Gaussian integrals do not 
result in additional prefactors. We recall that xi is understood to have a positive imaginary 
part, which is sufficient to ensure the convergence of (|i~8|). 

The ensemble average can now be performed immediately by substituting ([0]) into ( |l5|) and 
completing the square in the exponent. Interchange of the orders of the C- and ^-integration 
is justified because the C-integral is uniformly convergent in <j>. In order to proceed, we rewrite 
the terms appearing in the exponent of fllgj): 

1 k 

A • ct . <t> B = E(-!) eTr c] ■ ® 4 P ) . (2°) 

e=0p=l 

and the complex conjugate equation for 4> B ■ C • 4>a- Here, (f)® and ^ are commuting and 
anticommuting vectors of length N which represent the components of <f>A and (pB ■ The C- 
integral can now be performed using 

1 



d[C]exp -iVS 2 Tr &C - i Tr (C*X + CY) =iV- 1 exp 



NY? 



Tr(Xy) 



(21) 



where N^ 1 = J d[C] exp [-iVS 2 Tr CflC~\ . (For N f = we have N c = N D in @.) The matrices 
are arbitrary complex N x N matrices. Thus, the generating function ([l^) becomes 

1 fc 



• • • ,ik) 



d[<p A ]d[(f)B] exp 



Af£ 2 



ee'=0 



x exp 



E(-i) e (E( x -j 

p=i 



pp'=i 



lp + ^Bp^Bp) ~ ^ T {^Ap4>Bp + 4>Bp<l>Ap)) 



e=0 



(22) 



As a result of the ensemble average, we obtain a fourth-order term in the exponent which 
should be decoupled using a Hubbard-Stratonovitch transformation. In order to accomplish 
this, we rewrite the first exponential in (^) as exp (— Str (AB)/NY?) . The (k + k) x {k + k) 
super-matrix A is given by 



A 



lm 



b Al ■ <f>Am) 



(23) 



with a similar form for B. Now, we can 'undo' the transformation in (21) by introducing a 
non-hermitean complex (k + k) x (k + k) super-matrix, a. 

1 



exp 



NY, 2 



Str (AB) 



d[a] exp -A^S 2 Str (cr 1 " • a J - % Str (cr 1 " • A + cr • B 



(24) 



After writing Str (cr"!" - A) = <f>\- (c^ <E> In) ' 4>A and Str (cr- A) = (f> B ■ (a<S> ljv) ■ <Ab, and performing 
a shift of integration variables given by a— > cr + x — j, we can express the generating function 
for the correlation functions as 

Z k (j u ...,j k ) = J d[<j> A ]d[(t> B ] J d[a] exp{-AT£ 2 Str [(crt + x - j) • (cr + x - j)] } 

x exp {-i<f> A ■ (cr t <g> Iff) ■ 4> A - icj) B ■ (a ® 1 N ) ■ 4> B - mT{(j)\ ■ (f) B + <fi B ■ </>a)} • (25) 



Here we still use the definition ([17]) for x and j, however from this point on x p and j p in ([jj 
are not meant to be multiplied by the N x N identity matrix. 

At this point, we would like to change the order of the integrations in ( p5| ) and perform 
the Gaussian integrals in cj>. This would result in a superdeterminant involving the matrix 
a. However, this is possible only if the </>-integral is uniformly convergent in a, which is not 
the case if cr is an arbitrary complex matrix. However, it can be shown in general H3] that 
it is possible to choose a certain non-compact parametrization of the cr-variables that ensures 
uniform convergence. An explicit construction for one- and two-point correlation functions has 



been given in jig, 25]. The parametrization in the case of the one-point function will be discussed 
in great detail in the next section. 

Postponing further discussion of these points, we arrive at the following expression after a 



change in the order of the integrals in ( 25 ) 

Z k (j) = J d[a\ exp (-iVS 2 Str (o ] + x - j) • {a + x - j)) sdet" 



N 



a* ttT 
ttT a 



(26) 



This partition function, with N appearing only as an overall parameter, is amenable to a saddle 
point approximation. 

Before we continue with the evaluation of this partition function, we will show in the case 
of the one-point function that the imaginary part of the partition function can be obtained 
by replacing the non-compact parametrization of (^) by a compact parametrization. We will 



conjecture that all higher-order multipoint correlators Rf~, as opposed to the defined in (12), 
can also be obtained by using a compact parametrization with the cr-variables parameterized 
according to 



a = USV- 1 . (27) 

This choice of integration domain allows us to utilize the Itzykson-Zuber-like integrals developed 
earlier p4|, 35, 3£, 37 1 for the integration over the super-unitary matrices. Unfortunately, we 



have not been able to construct a rigorous proof of this statement. 

4 The one-point function 

In this section, we wish to make the point that the spectral density can be obtained from a 
compact parametrization of the cr-matrix rather than the non-compact parametrization which is 
required for uniform convergence, and which justifies the interchange of <p- and a-integrations. 
Since fermionic integrals are always finite, convergence problems arise only from the integrations 
over the boson-boson block of the a matrix. 

The uniform convergence of the (^-integrations can be achieved if we perform a Hubbard- 
Stratonovitch transformation using the identity 

„ „ 1 /"OO [-OO 

e~ a +b =— ds adae -« 2 -2iaacoshs-2ibasinhs ^ ^ 

7TZ J — oo J — oo 

where a is real positive and b 2 — a? has a negative real part. When expressed in terms of 
the real and imaginary parts of <jbb = 0*1 + i<?2, this identity corresponds to the non-compact 
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parametrization introduced in [25]. This parametrization is given as 

o\ = (a — ie) cosh s/E , 

02 = i{a — ie) sinh s/S , (29) 

where both cr and s run over the real line. Note that this parametrization covers only half of 
the (<7i,<T2) plane. (The parameter a introduced here should not be confused with the matrix 
a defined in the previous section.) We wish to contrast this parametrization with the compact 
parametrization (i.e., in polar coordinates), 

o\ = a cos 9 /Tj , 

£j 2 = a sin 9 /£ . (30) 

After performing the integration over the Grassmann variables in the cr-matrix, the partition 
function can be written as 

j\j~2 poo roo roo r2n 



J\j4 roo roo roo rzir 

Z(j) = — / ado / ds pdp d<pF((o - ie) 2 , p 2 ,t 2 ) 

J-oo J-oo JO JO 



-oo "'O ./0 

(2 2 \ N 

P + I | e -N[o 2 +p 2 +2(x- j)Y.(o-ie) cosh s+2ixY,pcos tp+Y, 2 {(x~jf -x 1 )} /^l) 

t 2 — (a — ie) 2 J 1 

where t = irTT,, p, and if parameterize the fermion-fermion block of the a matrix. We have 
displayed the ie dependence of (|3l]) explicitly. The convergence of the s-integral is guaranteed 
by choosing x — j to be pure imaginary. The definition of the function F((a — ie) 2 ,p 2 ,t 2 ) and 
additional details can be found in [25|. 

We perform the cr-integral first by saddle-point integration. The cr-integrand has poles at 
±t + ie, which allows us to deform the integration contour into the lower complex half-plane. 
The saddle-points in the a integration are given by a = ±i\/l — t 2 , and thus only the saddle 
point with the negative sign can be reached by deforming the integration contour. The resulting 
integral can be analytically continued to x — j just above the real axis. 

We now consider the imaginary part of the partition function. The ip integral gives rise to a 
Bessel function, Jo. Thus, the only contribution to the imaginary part comes from the factor 

/oo 
ds exp[-2A^(x - j)Ea cosh s] . (32) 
-oo 

For positive x, we change integration variables according to s — > s — iri/2. Then, the integration 
path can be deformed into the integration path I shown in Fig. 1, and S is given by 

S = J ds exp[2iiV(x — j)Ecj sinh s] . (33) 
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Im(s) 



Rc(s) 



II 

Figure 1: The integration contour for the imaginary part of the generating function. 

The imaginary part of S is given by half the sum (note the arrows) of the contributions along 
the integration paths I and II in Fig. 1. The integral along the parts of the integration contour 
parallel to the real axis cancel, and we are left with 



1 \ m 

Im S = —(S — S*) = Im - / dsexp[2iN(x — j)£<rsinhs] . 

Z Z J —ni 

If we change integration variables, s = i9, we find 



Im S 



d9 cos[2iiV(x - j)Ea sin 9] 



(34) 



(35) 



In effect, this implies that the imaginary part of the partition function can be obtained by replac- 
ing the (a, s) parametrization required for uniform convergence by a compact parametrization 
such as (^). With this parametrization, the partition function is given as 



jy2 poo r2ix 

Z c (j) = I a da I d9 



vr^ Jo 



pdp I d(pF((a - ief , p 2 ,t 2 ) 



p 2 + t 2 



N 



-N(a 2 +p' 2 +2(x-j)Y,(o-ie) cos 0+2jx£pcos tp+Y, 2 ((x-j) 2 -x 2 )) 



■ (36) 



\t 2 -{a-ie) 2 ) 

In this form, the imaginary part of the partition function is determined by the ie in the singu- 
larities of the pre-exponential factor. Therefore, the imaginary part is even in a. This allows 
us to extend the a-integration from — oo to oo with the introduction of a factor 1/2. Such an 
extension of the ^-integration to the full real axis is required if we are to perform this integral 
by saddle-point integration. This integration can be carried out in the same fashion as in the 
non-compact case and yields the same result for Im Z c (j) as obtained by substituting the final 
result for ImS (j35|) back into Im Z(j) defined in (|3ll), 

12 



By contrast, the real part of the integral is odd in a, and it is not possible to extend the 
cr-integration over the entire real axis. Consequently, the cr-integration for the real part of Z(j) 
cannot be performed by saddle-point integration if a compact parametrization is adopted; only 
the imaginary part of Z(j) can be obtained in this fashion. A similar phenomenon was observed 
in |^6| for both the one- and two-point functions. 

This leads us to the following conjecture: The spectral correlation functions, which are given 
by the products of the imaginary parts of the resolvents at different points, can be obtained by 



parameterizing the cr-integration in (26) as 



a = U(S- ie)V~ x , (37) 

where U and V are taken from the super-unitary group as discussed in the next section and 
where S is a diagonal matrix. This conjecture is central to the remainder of this paper. 

5 Integration over soft modes using Itzykson-Zuber integrals 

We now return to the evaluation of the generating function for correlations. One possible 
method for the evaluation of (|26|) would be to integrate out anticommuting degrees of freedom 
first. This is simple for a small number of Grassmann variables when only a small number of 
terms contribute to the integral. However, even for k as small as 2 this approach is impractical 
H3] , and for larger values of k it is virtually impossible to perform the integrals this way. 

Rather, we adopt a technique developed in [p2j ], where it was applied to a study of the 
Gaussian unitary ensemble. This technique consists essentially of separating the eigenvalue and 
angular coordinates of a by diagonalization and then integrating over the angular coordinates 
using a super-symmetric analogue of the Itzykson-Zuber integral. The remaining bosonic eigen- 
value integrals can then be carried out either exactly (e.g., for the case of zero temperature 
considered in section 5) or in saddle point approximation (e.g., for non-zero temperature as in 
section 6.) The main virtue of Guhr's technique is that it preserves the determinantal structure 
of the correlation functions, which makes it possible to obtain all correlation functions at the 
same time. 

We begin by reminding the reader of some relevant properties of super-unitary matrices. As 



in the case of an arbitrary complex matrix [47], an arbitrary super-matrix can be diagonalized 
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by two super-unitary matrices 



a 



U^SV 



(38) 



where S = diag(s5, . . . , sS; is\, . . . , is\) is a diagonal matrix which can be taken to have non- 
negative real entries and (U, V) parameterize the group (U(k\k) x U(k\k)) /[U{l)) k+k . (Here, 
U(k\k) denotes the super-unitary group) With this parametrization, the integration measure 
can be written as [] 

d[a] = d[S 2 ]dp(U, V)B 2 (S 2 ) , (39) 

where dfj,(U, V) is the super-invariant Haar measure, d[S 2 } = Ul=oI\Li d (sf) 2 - The Jacobian of 
the transformation is given by the square of the Berezinian |J6] . When the bosonic and fermionic 
blocks are of the same size (i.e. for Nf = 0) the Berezinian can be simplified to [j32f 



B(S 2 



det 



1 



(40) 



- p,m=l,...,k 

In terms of integration variables defined by the polar decomposition (|3^) (with an infinitesi- 
mal negative imaginary increment included in the eigenvalues) , the generating function ( |26| ) can 
be rewritten as 



ZfcCj) = / d[S 2 ] B2 (S 2 ) sdet" 

x j dii{U, V) exp (-iVS 2 Str (a + x - j) t (cr + x - j)) , 



S irT 
nT S 



where we have used the fact that 

sdet 



0"^ ttT 
ttT a 



sdet 



S TtT 
ttT S 



(41) 



(42) 



The group integral in (^) has precisely the form of an Itzykson-Zuber integral over complex 
super-unitary matrices. For arbitrary complex super-matrices a = U~ l SV and p = XJ ~^RV' , 
one finds by an application of the heat kernel method |34|, |35], [36], |37]] that 

^Vje-^IM'Ml = 1 det7( S °^)det 7 (^, r i) ^ 



(k\) 2 B(S 2 )B(R 2 ) 

where (s°; is},) and (r°; irV) (with p = 1, . . . , k) denote the eigenvalues of a and p, respectively. 
The quantity 7(5, r) is defined as 



7(5, r) 



exp 



2 2 
-s — r 



I {2sr) . 



(44) 



1 Note that our normalization of d[a] and dn(U, V) in (^3|) differs by a factor 2 2fc from ref. 
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For a derivation of this integral we refer the reader to [36]. In this reference contributions from 



Efetov-Wegner terms [ 33 1 are discussed as well. 



In (|4l|), the matrix p is diagonal with diagonal elements R given by 

R = diag(xi . . . ,x k + j k ;xi - ji, . . . ,x k - jk) , (45) 

and the matrix a is as given in (|38[). It is easy to verify that 

1 k 

Bfjpj = l[4x p j p x(l + 0(j)), (46) 

which enables us to carry out the differentiations with respect to the source terms. We observe 
that, as a result of (|40|), the integrand of the generating function ( |H|) factorizes into products 
and determinants of k x k matrices. By renaming the integration variables sp it can be shown 
that all k\ terms in the expansion of the determinants of j{Sp,r^) and 7(s„j,?"„) are equal. The 
only remaining determinant is the Berezinian (|40|). The products can then be absorbed into the 
determinant, which leads to the following determinantal structure for the correlation functions 

R k (xi,...,x k ) = det T,K N (x p E,x q E) , (47) 

L Ip,q=l,...,k 

where the dimensionless kernel i£j\r (£>??) is given by 

x e -N(r*+ s *+(e-v 2 )m Io (2NrC)I (2Nisrj) . (48) 

In this equation, we have introduced a dimensionless temperature 

t = ttTY, . (49) 

We have shown that all correlation functions follow from a single two-point kernel. This is the 
main virtue of the application of Guhr's supersymmetric method. 

The spectral correlation functions R(xi, ■ ■ ■ ,x k ) follow from the discontinuities of xi across 
the real axis. The imaginary part in (|48D arises as a result of the — ie term in the One can 
easily convince oneself that the spectral correlators preserve the determinantal structure (|47|). 
We thus find 

R k (xi, . . . ,x k ) = det[EK N (Ex p ,T,x q )] m=lt k , (50) 
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where 



8N 2 nr- f°° , f°° , rs t 9 ' 1 xJV 

Im 



8iV f°° f°° rs / \ j 

^v(£,r?) = v 7 ^? / dr ds^—^(s 2 + t 2 ) 

n Jo Jo r 2 + s 2 \ ) 



10 Jo 

-2 i '1 , /e2_„2\ 



-(r-ie) 2 + t 2 , 

x g-jv^+^+C^-VJ/a)^ (2iVr ^ /o (2iVi S7? ) . (51) 



Both in ( |48| ) and (^) the (x, y)-dependence has been chosen so that the kernel is symmetric in 
x and y when t = 0. 

Both at i = (section 6) and t ^ (section 7) the two integrals can be separated by the 
Feynman method, 

' N [°° dae- Na( - r2+s2 K (52) 



r 2 + s 2 Jo 



Most of the results leading to ([47p can be generalized immediately to an arbitrary number of 
flavors. The main modifications are that the fermionic blocks will be of size k + Nf rather than 



k, and that the supermatrices x and j in (17) will have Nf additional zero blocks. However, 
to the best of our knowledge, it is not possible to write the Berezinian as a determinant as in 
(flO|), and therefore the determinantal structure of ( p7|) is lost. This means that it is no longer 
possible to express the correlation functions in terms of a single kernel. 

On the other hand, work on both Wigner-Dyson random matrix models and chiral 
random matrix models with arbitrary unitary invariant potentials [^] using the orthogonal 
polynomial method shows that the correlation functions still possess a determinantal form. 
This suggests that the present approach might be modified to arbitrary Nf as well. 

6 Exact evaluation of correlation functions for T = 



In this section, we will show that kernel < p\\ ) at T = reduces to the usual Laguerre kernel 
obtained by means of the orthogonal polynomial method 

After the introduction of new integration variables by u = r 2 and v = s 2 , the kernel can be 
written as 

K N (x,y) = -—^ e - N ^-y 2 V 2 H dal u (x)l v (y) , 
vr Jo 
poo / 1 \ N 

IJx) = / dnlm J (2Nxy/u) e~ N( - 1+a)u , 

Jo \-u + ieJ 

roo 

T v (x) = / dvv N I (i2Ny^J) e ~ N ( l +^ . (53) 
Jo 
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In order to evaluate X u , we use the identity 



Im (u — ie) = ir 



(_l)jV-l d N-\ 

(N - 1)1 du N ' 



(54) 



in (53) and apply partial integrations in order to eliminate the derivatives of the delta function 
As a result, we find 



I u (x) 



(N-iy.du*- 1 



J (2Nxy/v) € 



-N(l+a)u 



(55) 



u=Q 



which precisely describes the derivatives of the generating function for the Laguerre polynomials 
(See eq. flO])-) We thus obtain 

vr(iV(l + a))^ 1 / ' 



2 u (x) 



NAT-l f Nx , 

-Ln-i 



(56) 



(N - 1)1 " X V 1 + a , 
where are the Laguerre polynomials. The integral j v is well-defined and can be found in the 
literature. (See eq. (|A.2[) .) The result is 

) / Nn 2 \ 

(57) 



e / Nv 



' (N(l + a)) N+1 ^- \l + a / 
After making a change of variables to z = (1 + a) -1 , we find for the two-point kernel that 

K N (x,y) =2N 2 ^e- N ^ 2 - y2 V 2 ^ dze~ N y 2z Ln^Nx 2 z)L N (Ny 2 z) 



To further simplify (|58|), we will need the identity (see eq. ( |A.5|) ) 

?7 — £ 

The integrand is now a total derivative, and we reproduce the well-known result 21] 



K N (x,y) = 2N^f 2 e 

x z — y z 



xy N(x 2 +v 2)/2 



L N ^(Nx 2 )L N (Ny 2 ) - L N {Nx 2 )L N ^{Ny 2 ) 



(58) 



(59) 



(60) 



This justifies our claim that the spectral correlation functions can be obtained by using a compact 
parametrization for the c-variables. 

It can be shown from the asymptotic properties of the Laguerre polynomials that spectral 
correlations in the bulk of the spectrum are given by the Gaussian unitary ensemble. The result 
for the microscopic region, x = Nx ~ 0(1), follows from the asymptotic form of the Laguerre 
polynomials, 



lim L r , 



n 



M^) , 



(61) 
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which can be used after rewriting (|30|) with the aid of recursion relations for the Laguerre 
polynomials. As a result, we find the microscopic kernel 

K s (x,y) = l\^^K N (^,^) = _v^(£j (2y)J 1 (2x) - yJ (2x)Ji(2y)) , (62) 

which agrees with results obtained previously |38|, |4(| . 
Finally, the microscopic spectral density is given by 

ps(x) = lim Ks(x,y) (63) 
y >x 

= x(J 2 (2x) + J 2 {2x)) , (64) 
which is also in complete agreement with previous results. 

7 Correlation functions at nonzero temperature 

In this section we evaluate the microscopic limit of the imaginary part of the two-point kernel 
at nonzero temperature and show that, up to a scale factor, it is in agreement with the zero 
temperature result in the limit N — ► oo. To this end, we perform the u and v integrals by a 
saddle point approximation. However, the saddle-point approximation to ([51]) suffers from the 
difficulty that u = — 1 + t 2 and v = 1 — t 2 . As a result, the pre-exponential factor diverges at 
the saddle-point. The aim of the transformations performed in the first part of this section is 
to eliminate this factor. 

We begin by separating the r and s integrals (^) according to the Feynman method (|52"|), 



The Feynman parameter, a, is replaced by the new integration variable (3 = + a). After 

rescaling r — > fir and s — * /3s, we find 

Kn{x,u) = — \fx~^j [ j3d[3 [ rdr [ sds ' ' ^ 



vr v y J " ^ Jo J y-( r -ie) 2 + tyf3 2 J 
x e -jv(r 2 + s 2 )+(£ 2 -s 2 )/2 Jq piprx) J {2(3sy) , 

(65) 

where we have also written the modified Bessel functions in terms of ordinary Bessel functions. 
Next, we express the product of the Bessel functions as a derivative of the microscopic kernel 
(|62^) according to the following remarkable identity (see eq. (|A.7| )) 



2/3^J (2/3x)J (2/3y) = A-(3K s (px, fy) . (66) 

dp 
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This identity can be derived from eq. ( |59|) using the asymptotic limit (|6l]) of the Laguerre 
polynomials. After insertion of this identity in (|65| ) and partial integration with respect to /?, 
we find 

K N (x,y) = -™1 f X dr r ' dsy/Fse-W+^-PU* 



vr Jo 



s 2 + t 2 \ N f 1 d ( s 2 + t 2 lj3 2 

r-ief+e ) K s(irx,8y)-J Q d(3K s (i(3rx, fit®- ( _ (f _ - ^ 



The second term in this equation can be simplified further. We differentiate with respect to 
s and undo the change of integration variables at the beginning of this section, i.e., r — * r/(3, 
s — ► s//3, and a = (1 — [3 2 )/[3 2 . Finally, we perform the integration with respect to a and 
obtain cancellation of the factor r 2 + s 2 which results from the differentiation. We find that the 
integrands of the two terms in (67) differ only by a factor —t 2 /(s 2 + t 2 )(—r 2 + t 2 ). Thus, this 



equation can be rewritten as 



oo 



K N (x,y) = -— l~ dr l~ ds^se-^^ 2 -^' 2 



vr Jo jo 



s 2 + t 2 \ N ( t 2 



(r-ie) 2 +t 2 \ (s 2 +t 2 )(-r 2 +t 2 ) 



Ks(irx, sy) . 

(68) 



The saddle point evaluation of K should be performed separately in the microscopic limit 
and in the bulk of the spectrum. In the latter case, the asymptotic forms of Bessel functions 
enter in the saddle point equations. However, this is not the case in the microscopic limit with 
x fixed in the thermodynamic limit. The saddle-point approximation in the microscopic limit is 
particularly simple. At the saddle point, we find 

r 2 = -1 + t 2 , 

s 2 = 1 - t 2 , (69) 

and both the second derivatives with respect to r and s are equal to 4iV(l — t 2 ). Note that 
the f 2 is outside the integration domain. As discussed in section 4, the imaginary part of the 
integrand of K^{x,y) is an even function of r, which allows us to extend the integration range 
from — oo to oo at the cost of a factor of 1/2. It is then clear that the integration path can be 



deformed to reach the saddle point —iyl — t 2 . We cannot extend the integration path for the 
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evaluation of the real part of K^(x,y), and it is not clear how the integrals can be performed 
by a saddle-point method. 
As a result, we find 

^^ N 44 ) = iCKs ^ ) - (7o) 

For convenience, we have introduced the scaling factor ( = y/l — t 2 , which gives the temperature 
dependence of the spectral density in the neighborhood of A = 0. We have shown that, up to 
this rescaling factor, the kernel Kn is independent of the temperature. This constitutes the 
central result of this paper. 

8 Conclusions 

In this paper, we have shown that all correlations of the eigenvalues near zero, measured in 
units of the average spacing, are independent of temperature deformations of the unitary chiral 
random matrix model. This result extends previous work on the microscopic spectral density. 
Together with other recent work on the universality of correlation functions with respect to 
deformations that preserve unitary invariance, this firmly establishes the universality of the 
complete eigenvalue distribution in the neighborhood of A = 0. 

It is our conjecture that the correlations of lattice QCD Dirac eigenvalues near zero virtuality 
are in the universality class of the chiral Gaussian Unitary Ensemble (chGUE). This conjecture 
has been supported by lattice simulations of the average microscopic spectral density via the 
valence quark mass dependence of the chiral condensate. In view of the present results, it would 
be interesting to study the correlations of lattice QCD Dirac eigenvalues in the neighborhood of 
A = 0. Our prediction is that such correlations are given by the chGUE. 

The present results were obtained by a generalization to the chGUE of the supersymmetric 
method developed by Guhr for the Gaussian Unitary Ensemble. The strength of this method is 
that it preserves the determinantal structure of the correlation functions. As usual in the super- 
symmetric formulation of random matrix theory, this method also requires a proper parametriza- 
tion of the integration variables. Following work by Wegner and Efetov, it was believed that 
hyperbolic symmetry was an essential ingredient for the parametrization of the integration man- 
ifold. The surprising feature of the present method is that all spectral correlation functions can 
be obtained from a compact integration manifold. One reason might be that, because of UaO-) 
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symmetry, the resolvent satisfies the relation 



G(x + ie) 



G(x - ie) 



(71) 



and that all spectral correlation functions can therefore be obtained from a generating function 
that does not involve infinitesimal increments of opposite signs. However, this does not explain 
why Guhr's method also allows for a compact integration manifold in the case of the GUE. 
Clearly, more work is needed to address this issue. 

Our results are based on the choice of a compact integration manifold. We have provided two 
important pieces of evidence supporting this choice. First, a detailed analysis of the imaginary 
part of the generating function of the one-point function shows that a non-compact integration 
domain can be transformed into a compact one. Second, a compact parametrization of the 
integration manifold reproduces the exact correlation functions of the chGUE. However, the 
ultimate justification of this change of integration variables remains an open problem. We hope 
to address this question in future work. 
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Appendix A: Some useful identities 

A generating function for the generalized Laguerre polynomials is given by 



The Laguerre polynomials are defined by L n (x) = Lfj °(x). A closely related integral is given 



The Laguerre polynomials satisfy the following remarkable identity 

nL n _ x {zrf)L n {zO = (n - l)L n _ 1 (zr ? )L n _ 2 (zO + e zr >^ [ze- zr >L n ^(zr])L n ^(zO} • (A.3) 
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by 





By a recursive application of this relation and the Christoffel-Darboux formula, 

n-l 

E L k(OL k (v) = 7 [L n ^(OL n ( V ) - L n (OL n ^( V )} , (A.4) 

we find 

e^ Ln{zv)Ln _ l{z0 = de~* (L n -M)L n (zy) - LnWL^jzr,)) 

dz 7] — 4 

From the asymptotic from of the Laguerre polynomials, 

lim n~ a L a n {-) = x- a ' 2 J Q (2v^) , (A.6) 



we obtain the following relation for Bessel functions 



2Py/xyJo(2Px)M2Pv) = ^(3K s (f3x, py) , (A.7) 



where the Bessel kernel K$ is defined as 



K s (x,y) = -4^-2 (xM2y)J 1 (2x)-yJo(2x)J 1 (2y)) . (Ai 
x L — y z 



This kernel can be obtained from the Laguerre kernel ( pCf ) with the help of the asymptotic result 
( |A.7| ) after rewriting the Laguerre polynomials in the same order by means of the recursion 
relation 

(a + n)L a n _ x = xL« +l - (x - n)L* . (A.9) 
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